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Abstract 

A direct algorithm is suggested for the computation of a linear output feedback for a multi 
input, multi output system such that the resultant closed-loop matrix has eigenvalues that include 
a specified set of eigenvalues. The algorithm uses deflation based on unitary similarity transfor- 
mations. Thus we hope the algorithm is numerically stable, however, this has not been proven as 
yet. 

1 Introduction 

Deflation is a technique that has been efficiently used in the solution of the standard eigenvalue 
problem of a matrix A as well as other eigenvalue related problems. According to this technique once 
an eigenpair of A is computed, we continue the process with a matrix that possesses only 

the remaining eigenvalues of A, and possibly a zero eigenvalue in the place of X\ . In this way we are 
left to solve a smaller problem. Deflation can be accomplished by a variety of algorithms. Some of 
them are, Hoteling’s deflation, Wielandt’s deflation, deflation based on similarity transformations, 
deflation by restriction, etc. An excellent review of the first three methods can be found in [7, pp. 
584-600], whereas deflation by restriction can be found in [5, p. 84]. Lately, deflation has been 
used in the solution of eigenassignment problems. For example, Wielandt’s deflation is used in [6] 
to solve the partial eigenvalue allocation problem with state feedback for continuous time systems. 
In this paper we will be concerned with deflation based on unitary similarity transformations, and 
how it can be used in the pole assignment problem with output feedback, or as it is often known, 
the eigenvalue allocation problem with output feedback (MEVAO). 

In section 2 we define the MEVAO. In section 3 we define transmission zeros and discuss how 
they affect the MEVAO. In section 4 we give an algorithm for the solution of the MEVAO. Finally 
in section 5 give some numerical examples to demonstrate the performance of our algorithm. 
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32.776 

3 

71.334 
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4 
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5 
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Table 11: Optimal weight using sequential appproxi matrons solution 
algorithm with a physical state-space realizat ion, for direct output feed 
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Lower case Greek letters represent scalars, upper case Roman represent matrices, while lower 
case Roman represent column vectors and indices. Superscripts T, H indicate transpose and con- 
jugate transpose respectively. The notation k = i(r)j means that k takes all the values starting 
from i until j with step r. e t represents the ith column of the identity matrix I, and A(-) the set of 
eigenvalues of a matrix, counting multiplicities. R(-), N(-) will represent the space spanned by the 
columns of a matrix, and the null space of a matrix respectively. .R,C will represent the set of real 
and complex numbers respectively. 

2 The Problem 

The MEVAO problem that we will attack in this paper may be defined as follows: 

We are given a controllable and observable system 


*(*) = 

Ax(t) + Bu(t) 

(1) 

y(t) = 

Cx(t) 

(2) 


where A G J2 nXn , B G i? nXm , C G i£ pXn , rank (B)=m, rank(C)=p, ar(<) G JR n is the state 
of the system at time t and u(t) the input at time t. We are also given a self conjugate set 
of eigenvalues A,, i = l(l)r with r=min{n,m+p-l}, we need to compute K G R mXp such that 
A(A — BKC ) D {At | i = l(l)r}. As a result, a linear output feedback u(t) = ~Ky(t) may be 
computed that will make (1) become 


x(t) = (A - BKC)x(t). 

The resulting system will have desirable properties for the contror engineer. 

3 Transmission Zeros 

In the remaining of the paper unless otherwise stated, we will assume without loss of generality 
that m < p. 

Transmission zeros ( trzs ) play an important role in the MEVAO, thus it is vital that they are 
well understood. Our experience with the literature on trzs has been rather disappointing. In our 
search for a definition we have found several; some of them even contracting one another. Thus we 
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decided to resort to the physical motivation of trzs, and from that to derive a sensible definition. 
Some observations that may give some physical motivation to the concept of a transmission zero 
(of a controllable-observable system in our case) were found in [2], p.41. From this we concluded 
that, physically, a transmission zero ( trz ) of a system is a specific frequency of the system for which 
there exists a nonzero input that will yield a zero output. This basically is definition 1 below. In 
the following we will present three definitions of trzs, and we will prove their equivalence. The 
reason we present these three defintions is because the first is directly derived from the physical 
motivation of trzs, the second is very useful in matrix computations, and the third is frequently 
encountered in the literature. 

Definition 1: A number < 6 C is a trz of (1),(2) if and only if there exists nonzero input u that 
can yield a zero output y = G(£)u = 0, where 


G(s ) = C{sl - A)~ X B , with s G C 


is the transfer function of (1),(2). D 

The number of trzs of (1),(2) can be at most n - max{m,p } (see for example [2], p.65). The set of 
all trzs of (1),(2) will be denoted by Zt T . Often it is helpful to use the following definition of a trz. 
Definition 2: C € C is a trz of (1),(2) if and only if 

^ =0 with x = ((I — A)~ 1 Bu . 

□ 


3u ^ 0 : 



Frequently we encounter in the literature the following definition of trzs . 

/*iia \ 

W<7 


Definition 3: Let = 







\ ° i 

be the Smith-McMillan (S-M) form of G(C), where L{ Q, R{0 are polynomial matrices and we have 
assumed that G(s) has normal rank m. Then C is a trz of (1) if and only if there exists i = 1 , ■ • • ,m 

such that £i(C) = 0. D 

In the S-M form of G(() the following properties hold 
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<i(C)l«a(0|-|«m(0 and tf m (Q | ^m-i(C) I • • * I V’i(C) • 


Because of the above properties we see that there cannot be trz £ that will also be a root of VVnCC)* 
Since, if it were then e m (C) — 0 and C would have been cancelled in €m(0 I ^m(0- Hence C would 
not be a trz , which is a contradiction. The roots of V>»(C)> 1 = l(l)m are the poles of (1),(2), 
counting multiplicities. In our attempt to prove the equivalence of the above definitions we will 
be using the expression G(Qu with £ £ Z tr • One then may wonder about the existance of G(C)u 
when C is also a pole of (1),(2). The following lemma will help us clarify this point. 

Lemma 1: Let ( G C be a trz as well as a pole of (1),(2), then there exists u / 0 such that 
lim G{s)u exists. 

Proof: Let £ be a pole as well as a trz of (1), then from the S-M form of G(£) we have that 

3 (i,j) with i < j : ifc(C) = ej(C) = 

Let us now for illustration assume p = m = 3 and ^i(C) = 62(C) = 0* Then by taking u T = 
(0 ,r?,0) / Owe have the following from the S-M form of G(s). 

lim L“ 1 (s)G(s)Ji” 1 (s)tt = 


/ lim V W 


^l(s) 


V 

/ o 
o 


lim 


\ 




lim5& 


\ 

(o\ 


V 

J 

\ 6 J 


which is defined since ^(C) / 0. Since L(s), R(s) are nonsingular (they have determinants inde- 
pendent of s, see for example [7], p.21), limZr 1 ^), Y\mR~ l (s) are defined. 

s-+c 

Hence if we take u = lim R~ l (s)u then lim G(s)u is defined for every £ £ Z tr that is also a pole of 
® *■ C s— C 


( 1 ),( 2 ). 


□ 
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In the following when we write G(Qu with ( E Z% r also being a pole of (1), we will mean 

limG(s)u. Similarly when we write x = (£/ — A)~ 1 Bu we will mean x — lim(s/ — A)~ l Bu y when 
s— ►C *“ t, C 

C G Ztr is also a pole of (1),(2). 

Theorem 1: Definitions 1,2 and 3 are equivalent. 

Proof: To prove that definition 1 is equivalent to definition 2 we proceed as follows: 


C G Z iT {3 u / 0 


{3u/0 


]«/0 


^ C 6 4 


C((I — A)~ 1 Bu = 0} 

Cx = 0 with x = (£1 — A)~* Bu} 
A -(I B W x 
C 0 I l « 


= 0 with x = (£7 — A) 1 Bu 


To prove the equivalence of definitions 1 and 3 we proceed as follows: 


(eZ tr {3«/0 : G(C)« = 0} 


diag id£i 


3u/0 : L(() «= R (0 U _ o 


A /(0 


( 3 ) 


Since L(s), R(s ) are nonsingular for Vs € C, M(() must not have full column rank in (2). Hence 


(3) {3 * 6 {l,2,- -,m} 


def 3 
< \ * 


C G Ztr • 


€i(0 = 0} 


□ 


In the last theorem we silently assumed that the normal rank q of G(s) is m, or equivalently 

A -si B 


the normal rank n + q of P(s) = 


is n + m. However, theorem 1 holds even for 


0 


the degenerate case where q < min{m,p} (= m according to our assumption). To show this we 
proceed as follows: 
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First we prove Z iT — C using definition 1. Since q < m the following is true 


{VC € C 3u ^ 0 : G(()u = 0} © Z tr = C . 


( 4 ) 


To prove the equivalence of definitions 1 and 2 we see 

/ 


ill 


A -(I B 
C 0 


Z tr = C |VC G C 3u t 0 : 

To prove the equivalence of definitions 1 and 3 we see 


= 0 with x = (Cl — A) x 5u} 


Zt T — C { 


VC 6 C 3u ^ 0 : L(C) f 1 R(C)u = 0 


{VC € C 3* €{!,•••, m} : c,(C) = 0}. 


Hence another definition of trzs that is often found in literature, and according to which C e z tr 
if and only if £?(£) and P(() go below their normal rank, does not agree in the degenerate case 
with definitions 1, 2 and 3. This is so because not G C will make G(C) and P(() go below their 
normal rank, hence Z tz ^ C , according to this definition. In what follows we will assume that 
trzs are defined by definitions 1, 2 or 3. 


Definition 4: A number £ E C is a strong transmission zero ( strz ) of (1),(2) if and only if 
G(() = or equivalently €\(C) = 0 (in the S — M form of G(£)), or equivalently 


Vu#0, 

( A- Cl B \ 

H 

= 0 , with i = (C I - A)~ l Bu 

l 

C °) 

\ u ) 

4 


Theorem 2: ( 6 C cannot be assigned by output feedback if and only if £ is a strong transmission 
zero. 

Proof: Let G(Q / 0, then if ( 6 A(A) we ta k e K = 0 and £ is placed. If C ^ A(A) then there 
exists u 0 such that 


V = G(0 U ? 0 




x = ((I- A)- x Bu 
y = Cx/0 
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If we take K = — then we see that 


[</ -(A- BKC )] x 


(Cl - A)x + BKCx 


Bu-B 



= 0 


Thus the eigenpair (£, x) has been placed, which proves the necessity of the propositionn given by 
the theorem. To prove its sufficiency assume G(Q = 0, but nevertheless £ has been placed. Then 
there exists x / 0 such that for some K 


(A - CI)x - BKCx = 0 =► x-(A-CI)~ l BKCx = Q 


=> Cx-C(A-(I)~ l BKCx = 0 
=» Cx + G(QKCx = 0 
=> Cx = 0 . 


Hence from (A - £J)x - BKCx = 0 =>■ (A- (I)x - 0 => £ € A(A). This is a contradiction, since 
if £ (E A(A) the G(£) is not defined, however it is assumed that G (( ) = 0. □ 


4 The Algorithm 


We give an algorithm based on deflation by unitary similarity transformations. 

Step 1: Compute a feasible eigenvector x\ of A\ — B\K\C\ = A - BKC corresponding to, say 
Aj with || x\ [( 2 = 1. We will show in the sequel how we may compute such an x\. 

Step 2: Allocate the eigenpair (Ai,xi) to A\ — B\K\C\ as follows: Let V\ be a unitary transfor- 

( Ri \ 

mation such that V^CiXi = b\e i and B\ = {Uq,Ui) I a QR decomposition of B\, then 


(A\ — B\K\C\)x\ = Aixi 






Axzx-BxiKiVxyyfCxxx) = Aix x 

= (Ai - Ai/)xi , with K\ — K\V\ 
B\k\Si = (A x - A x /)x x , with ki = k\e\ 

[ * ) M, = ( VS{M ~ 

\ o ) ^ U( f (A 1 - Aj/Jx! 


( 5 ) 
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From (5) we see that if x\ was computed in step 1 so that x\ E ( A\ — Aj/)], and k\ is computed 

by solving 

* 1*1 = st 1 u?(a 1 -\ 1 i)z 1 w 

then (5) is satisfied and the eigenpair (Ai,xi) is allocated. An eigenvector X\ that satisfies 
x\ E N[f/ 1 r (Ai - Ai /)] will be called feasible. 

Step 3: Compute unitary <5i = ( X \,Q\ )• Hence from Q\e\ = => Q\X \ = Ci we see that Q\ can 

be a Householder transformation or a sequence of rotations. 

Step 4: Perform the unitary similarity transformation 


( x H \ 

! [A 1 -B 1 K 1 (V 1 h C 1 )](x 1 ,Q 1 ) 

\ Qi ) 

j [(A x x x - B x k x 6 x ) , (A X Q X - B X K X {V”C X )Q X )) (7) 



From step 3, fci has been computed so that A\X\ — B\ki6\ = AiXj. If we now set K\ — (A:i, # 2 ) 
and 


V X H C X Q X = 



then 



where A 2 = Q X A X Q X - Q? B x k x c% , B 2 = Q? B 2 . 

Step 5: Continue the allocations with A 2 - B 2 K 2 C 2 . □ 

There are two points that need to be clarified for the allocation of Ai in the above algorithm. One 
is, how to identify whether Ai is a strz or not, and the other is the computation of * 1 . 
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Lemma 2: If {Ax is a strz of (1) } => 6 X = 0. 
Proof: Let A x be a strz of (I) then 


V«x / 0 3x\ 



= 0 


=> C\X\ =0 => <5x = 0 


We may note however that the above lemma on its own is not adequate to identify strz . In the 
sequel we will compute our feasible eigenvector x\ in such a way that if 6i = 0 then Ax will be a 
strz . As we will show, computing xx in this way will also improve the numerical stability of our 
algorithm. To accomplish the above we proceed as follows: 

We know that 

i 1 = 11 C\X\ || 2 =* {<$x = 0 xx E N(Cx)} . 

Thus if we compute xx so that 


X! G N [U»(A l - X l I)) - N(Cj) 

when 

then we may safely conclude that £x = 0 => {Ax is a strz }. If the dimension of 

N [V?{A X - A,/)] - N(C,) 

is greater than one, we have some freedom in choosing our Xx. To see how we may exploit this 
freedom, we observe from (6) that if | #x | is unreasonably small then we may unnecessarily lose 
accuracy when we solve (6). To avoid this, let 

c? = (H 0 ,H i) | J 1 j 

be a QR decomposition of Cf, then R(//o) = R {Cf) and R(//j) = N(Ci). Assuming now 
II *i || 2 = 1 we get 


I 1 = 11 ||, =|| Zi HqXi ||,<|| Z, IUII H?x x ||, . 

But || HqX i |U = 0max(H qXi) = cosO , <T max (-) is the maximum singular value of a matrix, and 
6 is the smallest angle between R (Cj) and R(xi) (see[l]). Hence from 
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I h I^H z || 2 cos9 


in order to make | 6\ | as large as possible we need to make 0 as small as possible. Thus we will 
choose ari E N[f// / (,4i — X\I)] that forms an angle 9 with R(7/o) as close to zero as possible. The 
following theorem taken from [1, p.582] and modified to meet our requirements, will give us a way 
to compute 0 as it was required above. 

Theorem 3: Let Ho and W\ form unitary bases of two subspaces. Let also 

H^W X = YXX h 

be the singular value decomposition of HqWx , then the smallest angle between R(//o and R(W r 1 ) 
takes place between vectors HoYei and W\Xe \ . 

Proof: See [l] 

Theorem 3 suggests the following algorithm for the computation of x\\ 

Compute the QR decomposition of C^, B\, and [Ui{A\ — X\I)] H = (Wq, W\) 

Compute the singular value decomposition of V^W\ = YYiX h . 

Take x\ = W\Xe 

Finally to eliminate any doubt that 61 = 0 when Aj is a strz we prove the following lemma. 
Lemma 3: - Ai/)] C N(Ci) {Ai is a strz of (1),(2)} 

Proof: First we prove 



N [Uf{Ai - A,/)] = {(A x / - A 1 ) -1 fi 1 ux | «x G R m } . 
To do so, we see that, given u\ ^ 0 and computing x\ to satisfy 


( 8 ) 


(Ai — X\I)x\ = — B\U\ <==> 


( u" 


(A x - Ax I)x\ = 


-R 1 


[u f 

=> xx 6 N[t/"(Ax - Ax/)]. 


W 1 


Thus 


{V«1 € R m , (A 1 / - Ax) _1 5x«i G N[t/"(A, - Ax /)]} 
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{(Ax/ - A \ ) — 1 xi\ I «1 e R m ) C N[t/, H 04i - A ,/)] . 


( 9 ) 


Let X\ : U\{A\ - A \I)x\ = 0, then, since we can always compute 


u G R m : Ri^i — Uq{A\ — A il)x\ , 


we have 


Uff 

V? 


(y4i - AiJ)*x = 


from which we get 


-R i 

o 


ux => ix=(Ax/— A\) 1 B\u i, 


N[^"(Ax - Ax/)] C {(Ax/ - ^i)- 1 .Bitii | «x € R m } . (10) 

From (9), (10) we may derive (8). 

Suppose now that 


{Ax is a strz of (1),(2)} <=► {Vux € R m C\(\\I — A\) 1 5x«i — 0} 

<=> CxN[t// / (Ax - Ax/)] = 0 

N[lf^ (i4x - Ax/)] C N(C,) . 


□ 

The algorithm that has been described so fax in this section, can be used to allocate only one 
eigenvalue. We will show however that we may use it to allocate m»n{n, m 4- p — 1} eigenvalues. To 
do so we need to observe that only the first column of the current Ki is needed for the allocation 
of one eigenvalue. We also need to consider the fact that A (A, — BiKiCi ) = A (Af — Cj Kj Bj). 
Given these two points then we may use the following algorithm, which for illustration we describe 
for m = 2 and p = 3, thus 

/ x x x ^ 

K x = 


XXX 

XXX 


We allocate eigenvalues until the number of rows of the current Ki become greater than the number 
of columns. In our example this happens after the allocation of two eigenvalues, hence 


K 3 = 



150 


At this point we continue working with the transposed system — Cj Kj Bj . Since Kj has 
two columns we are able to allocate two more eigenvalues instead of just one that we would have 
allocated if we had continued with K 3 . In the general case we keep transposing until we run out 
of eigenvalues or columns. Thus the total number of eigenvalues that we manage to allocate using 
this algorithm is 


m + p — 1 = 4 . 

Note that by following this algorithm we also satisfy the condition m < p at the beginning of each 
allocation. This condition has been assumed throughout this paper. Note that m,p are associated 
with the columns of the matrix on the left of Ki or Kj and the rows of the matrix on the right of 
Ki or Kj respectively, rather than with B{ and C,\ 


5 Numerical Examples 


In this section we give two numerical examples to demonstrate the performance of our algorithm. 
The computation was performed on double precision (56-bit mantissa) using PC-MATLAB on a 
Toshiba T5100 which uses an 80387 coprocessor equipped with the IEEE floating point standard 
of arithmetic. In the computation below we show, up to 12 decimal digits of accuracy. 


Example 1: 


A = 


f 0.581314086914 
0.1667175292% 
0.353500366210 
0.836242675781 
^ 0.244094848632 


0.504058837890 

0.157394409179 

0.441650390625 

0.200820922851 

0.936126708984 


0.559921264648 

0.665222167968 

0.373367309570 

0.072296142578 

0.607955932617 


0.741333007812 

0.933609008789 

0.771942138671 

0.598373413085 

0.154357910156 


0.303421020507 ^ 
0.144439697265 
0.078048706054 
0.638671875000 
0.154678344726 1 


B t 


' 0.549163818359 
0.942672729492 
^ 0.613098144531 


0.843856811523 

0.008666992187 

0.065872192382 


0.178649902343 

0.239028930664 

0.300445556640 


0.409255981445 

0.142791748046 

0.576828002929 


0.595489501953 ^ 
0.671371459960 
0.726318359375 j 
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c = 


' 0.561660766601 
0.427825927734 
^ 0.455917358398 


0.332702636718 

0.648269653320 

0.134307861328 


0.355514526367 
0.666625976562 
0.49757385253 9 


0.279266357421 

0.572494506835 

0.212463378906 


0.531448364257 ' 
0.972076416015 
0.653732299804 j 


The eigenvalues to be allocated are 


{0.704589843750, 0.421768188476, 0.572113037109, 0.396102905273, 0.127380371093}. 


The K computed by our algorithm was found to be 


K = 


f -9.415631257941 
-36.155174309311 
^ 33.590525809043 


-1.520241736469 

-4.515055559417 

4.071903303950 


11.520696978251 

41.184483490595 

-36.923138190434 


The eigenvalues of A - BKC were computed and they were A(A — BKC) — 


\ 

/ 

{0.704589843749, 


0.421768188479, 0.572113037109, 0.396102905270, 0.127380371093}. 


Example 2: In the above example we had n = m + p- 1 , thus all n eigenvalues were allocated. 
However, if Tti p — 1 < n then n — (m 4- p — 1) eigenvalues of A — BKC will take values that our 
algorithm has absolutely no control over. The following example demonstrates this point. 


/ 

A = 


0.356336616284 

-0.210549376245 

-0.355939324751 


-0.356626507847 

2.165165719183 

-0.506195941988 


-1.198870998736 ^ 
-0.882324378697 
0.721098719251 y 


/ 

B = 

\ 


1.606955436561 

0.062823512419 

-1.611627967397 


-0.407338058168 ^ 
-0.595038063525 
0.616657720777 y 


C = -1.182820557312 0.343687857622 -0.357421120340 J 

The eigenvalues to be allocated are {-0.406648486670, -0.366406484747, 0.853280016421}. 
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The K computed by our algorithm was found to be 

( -2.053475126112 
-6.466096811958 

The eigenvalues of A — BKC were computed and they were A(A — BKC ) = {-0.406648486670, 
-0.366406484747, 1 .707616832510} . 

Our algorithm allocates one eigevalue at a time, thus complex eigenvalues need complex arith- 
metic. As a result, K may be complex. Investigation is under way to derive an algorithm that 
will allocate two eigenvalues at a time in a double step. In this way we will allocate a complex 
conjugate pair of eigenvalues in one double step using only real arithmetic. Hence K will be real. 


6 Conclusion 

We presented an algorithm for the pole assignment problem for multi-input, multi-output systems 
using output feedback. The algorithm uses deflation based on unitary similarity transformations 
and it allocates rmn{n,m + p — 1} eigenvalues. The same kind of deflation has been used in [3] to 
solve the corresponding pole assignment problem using state feedback. Since the algorithm in [3] 
has been proven to be numerically stable we hope the algorithm in this paper has the same property 
too. However, this needs to be proven by doing a rounding error analysis of the algorithm, and we 
plan to do this in the near future. 
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